Experimental identification of individual insect visual tracking delays in free flight and their effects on visual swarm patterns

Insects are model systems for swarming robotic agents, yet engineered descriptions do not fully explain the mechanisms by which they provide onboard sensing and feedback to support such motions; in particular, the exact value and population distribution of visuomotor processing delays are not yet quantified, nor the effect of such delays on a visually-interconnected swarm. This study measures untethered insects performing a solo in-flight visual tracking task and applies system identification techniques to build an experimentally-consistent model of the visual tracking behaviors, and then integrates the measured experimental delay and its variation into a visually interconnected swarm model to develop theoretical and simulated solutions and stability limits. The experimental techniques include the development of a moving visual stimulus and real-time multi camera based tracking system called VISIONS (Visual Input System Identification from Outputs of Naturalistic Swarms) providing the capability to recognize and simultaneously track both a visual stimulus (input) and an insect at a frame rate of 60-120 Hz. A frequency domain analysis of honeybee tracking trajectories is conducted via fast Fourier and Chirp Z transforms, identifying a coherent linear region and its model structure. The model output is compared in time and frequency domain simulations. The experimentally measured delays are then related to probability density functions, and both the measured delays and their distribution are incorporated as inter-agent interaction delays in a second order swarming dynamics model. Linear stability and bifurcation analysis on the long range asymptotic behavior is used to identify delay distributions leading to a family of solutions with stable and unstable swarm center of mass (barycenter) locations. Numerical simulations are used to verify these results with both continuous and measured distributions. The results of this experiment quantify a model structure and temporal lag (transport delay) in the closed loop dynamics, and show that this delay varies across 50 individuals from 5-110ms, with an average delay of 22ms and a standard deviation of 40ms. When analyzed within the swarm model, the measured delays support a diversity of solutions and indicate an unstable barycenter.


Introduction
Insects are model systems for resource-constrained micro air vehicles operating in group and swarm applications. Although these naturalistic swarms rely on limited sensory and neural feedback structures and lack a traditional engineered communication network, they achieve coordinated flight maneuvers in near proximity to neighbors and amidst changing neighbors in unstructured environments [1][2][3]. Vision is one of the few sensing modalities whose quantified bandwidth, range, and sensitivity could provide realtime feedback to modulate these flight paths, and the relatively large fraction of insect neural material dedicated to visual processing suggests that visual control may be an important tool for implicit communication [4][5][6]. The exact mechanisms that are used to facilitate this interconnection are not yet known. Many theoretical swarm models and analyses remain disconnected from experimental studies on naturalistic swarms, which limits the ability of these parallel efforts to inform each other [7,8]. In particular, the effect of latencies such as sensory processing delay in biological (or computational in robotic) swarm experiments has not been adequately connected to theoretical developments.
This study develops an experimental tool based on a flight arena equipped with a visual stimulus system (Stimulus design) and a real-time insect tracking routine called Visual Input System Identification from Outputs of Naturalistic Swarms (VISIONS) for insects in free flight (Tracking system), and uses this tool to experimentally identify the visual processing and feedback delay across different insects. The identified delays are measured at the individual insect level, and probability distributions are experimentally quantified for the measured delays ((System identification). The swarm level effect of these delay distributions are then connected to recent progress in theoretical swarm communications analysis and a implemented in a simulated swarm environment (Analysis & simulation). The analysis and simulation indicate that three characteristic patterns may be generated by a collection of insects visually navigating relative to each other with such delays, and identify the Gaussian delay distribution range needed to generate stable swarm (Results and discussion).

Experiments for insects' visuomotor response
The effects of environmental stimuli examining insect flight behavior and motions during visually-dominated behaviors like obstacle avoidance, landing on a wall or proboscis, and flower tracking have previously focused most on the role of ambient and external illumination levels [9][10][11]. By high-speed videography, it was observed free-flying fruit flies generating quick turns through modest wing alterations and discussed methods for a one-parameter open loop paradigm for free-flying insects under a visual display [12]. To adapt for rapid changes in body orientation, dipterans and bees adjust the amplitude and angle of their wing stroke asymmetrically [13][14][15]. Nocturnal sweat bees do not change their flight speed while landing in a tunnel [16]. Later work showed that diurnal insects like bumblebees and hornets reduce their flight speed and take progressively more circuitous paths as the light intensity decreases [17], however nocturnal bees Megalopta were able to maintain flight speed with decreasing light intensity, including a spatial summation mechanism to explain their nighttime heightened sensitivity [18]. Flower tracking investigations on Manduca sexta utilizing a moving flower at various frequencies indicate tracking performance during feeding and energy expenditure, revealing that flower movement direction has a stronger effect on tracking performance than the target's movement frequency [19]. When flying in the wind, the unsteady flight of a hawkmoth in the wake of a 3D printed robotic flower displays larger tracking overshoot and a reduced order dynamic system [20] and a flower tracking experiment was used to quantify the change in their flight behavior under various light conditions, with flower tracking behavior represented by a simple temporal delay at various light intensities [21]. By adjusting light intensity, a system identification approach was utilized to find a brightness-dependent delay term in the transfer function, resulting in a dynamic model for each Hawkmoth variant that included a combination of species-dependent scaling parameters and processing delays [22]. There is evidence to suggest the internal delay time may have a dependency on swarm sizes in mosquitoes tracking a moving marker [23] and by task; in robber fly predation and obstacle avoidance tasks, different time delays (30 and 90 ms) for proportional navigation and obstacle avoidance rule were observed [24].
While experimental research is beginning to understand the need to quantify internal delays due to insect sensing and feedback, these previous studies are limited to reporting a single average delay across all animals and do not yet account for the heterogeneity of delay across the population or the effect of such delays on neighbor-coordinated behaviors.

Camera based tracking technologies
Animal movement study has benefited greatly from availability of electronic cameras to reconstruct animal positions. Early automated tools for tracking insects like "Ethovision" and "Trackit" relied on color cameras and chrominance differences to segment foreground and background [25,26]. "Flydra" achieved real-time tracking and allowed the use of monochrome cameras without the need for chrominance-based segmentation by assuming a fixed background scene and using intensity based segmentation [27]. This technology began to move to other schooling animals, with parallel improvements in outdoor starling tracking applying trifocal cameras began to bring this technology outside the lab, assuming homogeneous sky through a fixed background assumption [28]. Other animals began to be tracked, like fish and mice position tracking from post-processed videos with "IdTracker", [29], and automotive traffic [30], both relying on a static background assumption. "Toxtrac" tracked salmon, zebrafish, and cockroaches using a fixed background assumption, and similarly relied on a static background with a well-lit subject [31], as did infrared tracking approaches on ants at 2 Hz update rate [32]. Further, visual tracking conducted in a naturalistic settings often results in lighting changes that negatively affect tracking approaches [29,31].
For studies involving an insect or collection of insects tracking a visual stimulus, the appearance of visual stimulus in the background is necessary for precise timing measurements (as will be discussed in the next subsection). Recent works allow limited noisy backgrounds [33,34], and many tracking problems are still completed manually [35]. Animal flight studies have progressed with the help of high-speed imaging systems that allow researchers to monitor the detailed kinematic behavior of animals during flight [36]. Using artificial markers on the wings wing deformations can be measured and quantified with high speed videography [37,38]. Because of the small wing sizes, rapid frequency flapping, and unpredictable movements in near proximity to the insect body, wing kinematics characterization and reconstruction remain a challenging task even with high-speed videography [39,40]. Recent improvements in high speed visual tracking digitizing wing motions have integrated mechanisms to handle dynamic backgrounds and occlusions [41] such as the moving average [42] and Kalman filter [43,44] techniques, and these approaches have each shown value in real-time visual tracking for robust visual tracking system identification.

System identification
System identification techniques are applicable to experimental manipulations in which an input and output stimulus pair must be used to discover a set of unknown underlying dynamics. These experiments generally involve prescribing inputs that perturb a system (sometimes referred to as a 'black box') with known input signals and recording the resulting output variables [45]. This identification problem is often solved in both time and frequency domains [46]. Frequency domain analysis includes the ability to generate non-parametric identifications, applicability to classical control-system design methodologies and modeling of flying qualities, noise robustness, and reduced dimensionality for model parameter estimations [47,48]. For a linear dynamic model estimate, a diverse set of methods based on frequency domain identification are available, such as equation error [49], filter and output error [46], and spectral estimation [50]. The relative maturity of such models has resulted in applications to insect flight behaviors such as chasing mates [51], saccade turns [52], obstacle avoidance [53], speed control [10], navigation [54], and flight stabilization [55][56][57].
Theoretical swarm models with delay Swarm models. Mathematical descriptions of collective motions of multi-agent biological systems, including bacterial colonies, slime molds, locusts, and fish, have become available in the last few decades [58][59][60]. These models range from continuous approximations and kinetic theory models to models at the individual level, with interaction mechanisms generally consisting of attractive and repulsive forces [61][62][63] that move the particles. Aerial insect sensory systems are implemented at the individual level, which lends itself to treating each biological or mechanical individual as a discrete particle. This structure is consistent with individual models like the evolvable heading Vicsek [64] and evolvable-velocity Cucker-Smale [65] models. Both of which are composed of a set of differential equations (each representing a single agent's dynamics) networked together by some idealized sensing and interaction rule. Rigorous proofs are available to show sufficiency of particular rules to lead to flocking-like behaviors showing velocity convergence.
Bifurcation analysis. In differential equation analysis, a bifurcation parameter can be used to understand when a dynamic system exhibits distinct changes in behavior, such as stability or topological structures such as saddle points [66]. Modulating a single bifurcation parameter may achieve distinct behavioral characteristics using pitch fork, Hopf, and transcritical bifurcation types [67]. Different swarming patterns of spacecrafts by attraction and repulsive potentials through dynamic system theory can be obtained [68,69].
Recent progress on theoretical delays. Despite the availability of rigorous proofs for sufficiency (e.g. [65]), the theoretical models have often relied on a high level of instantaneous connectivity that may be impractical in nature or robotic implementations. The theoretical effects of delay in such a network of interaction rules has previously been investigated primarily numerically by distributing delays among agents and applying a mean field analysis [70]. Only recently are these effects beginning to be understood theoretically by extending the Cucker-Smale model to include the effects of time delay, which provides a path to understanding the effects of sensor and communication delays in swarm networks. Communication weights and a Lyapunov approach are used to determine an upper bound on time delays for velocity convergence for the Cucker Smale model [71].
An upper bound on the size of time delays is now available as a sufficient condition for flocking behaviors [72]. These effects of fixed time delays have recently been applied to both first and second order swarm models with fixed time delays [73] in theoretical and numerical simulation. Bifurcation analysis applied to the swarm network problem with delay as the bifurcation parameter can yield insights into the delay structures that affect collective motion behaviors [74]. To overcome the existing mean field theory's predictability limitations, stability analysis of ring and rotating patterns in a delay coupled swarm is carried out utilizing the bifurcation parameter [75].

Contribution of this paper
Despite the progress in different tracking experiments on animals, there is a lack of archival literature directly connecting experimentally measured delays in swarming-capable animals to the theoretical structures available. It remains unknown to what degree visuomotor tracking delays vary across individuals in a population (or even whether they might be better be modeled as fixed or heterogeneous), and whether these individual delays support or conflict with current understanding of visual swarming. A primary contribution of this work is to address these gaps between theory and experiment by providing: (a) an experimental quantification of how individual insect visuomotor delays vary over population, and (b) a theoretical and simulation analysis of the effect of these delays on visually connected swarms, illustrating the achievable shapes and stability regions and the need for delay modeling.

Experimental data collection
We prepared a rectangular acrylic enclosure of 15 × 10 × 10 inch dimensions, as illustrated in (Fig 1a and 1b). Honeybees were collected from a research beehive and placed inside the enclosure at 72˚F. The experiment required constructing two pieces of equipment: a visual stimulus input and a real-time tracking system. The experiment moves the stimulus in the X axis of the world coordinate, while the insect flying in the working volume tracks it and we record its trajectories by the tracking system. The insects track the light intermittently rather than continuously (S1 Video), and this visual tracking study analyses the periods of tracking. Tracking sections S in a trajectory were identified automatically as follows, S ¼ ½DTðqÞ:::DTðq þ gÞ�; g � 30; q ¼ 1; ::: where k is the total number of points. ΔT(m) is the absolute difference of every pair of points defined by ΔT(m) = |T 1 (m) − T 2 (m)| < η σ . T 1 (m) and T 2 (m) are continuous trajectory point m of X coordinates of the stimulus and bee. The tracking tolerance η γ was used to specify the desired tracking accuracy. S ¼ ½S 1 1 ; S 1 2 ; :::S 1 g S 2 1 ; S 2 2 ; :::S 2 g � is the identified section containing g tracking point pairs where S 1 ð1:::gÞ and S 2 ð1:::gÞ are stimulus and bee points respectively.

Stimulus design
The visual light stimulus is provided via a two-sided rectangular LED display, each of them contains 16 by 32 (width × height) LEDs with a refresh rate of 100 Hz. A Linux computer and a Raspberry Pi microcontroller are used to generate a programmable image on the LED panel. Fig 1c depicts the two Adafruit dotstar LED displays, one 74aN174 level shifter, and a Raspberry Pi microprocessor, as well as their wiring. Our stimulus was a synchronized vertical green light with a height of 2 by 8 inches (width×height), and a brightness of 2200 lux (measured at the center point between the two opposing lights), as seen in Fig 1b. The stimulus is moved by a sum of sinusoidal frequencies that changed the location of the vertical bar along the world coordinate's X axis.

Tracking system
Three synchronized cameras (Flea mono cameras) record the bees' flight paths at 50-120 frames per second. Fig 1a shows the cameras set on three tripods with angular separations varying from 50 to 90 degrees. The raw images (1280x1024) collected by these cameras are transported to a central Linux computer through USB3 connections. The ROS platform's Python-based program runs in real time and may also be used for recording and analysis offline. The main flowchart for the tracking system is shown in Fig 2, and the major processes are discussed below. 2d centroids detection. After obtaining raw photos from the cameras, we distinguish the foreground from the background. When the background varies slightly, the moving average technique is appropriate. The initial background is found by taking average of first few frames. Then, the background E t (i, j) for each new frame at time t is updated as where i and j denote the pixel coordinates that change in every frame, and F t (i, j) is the current frame at time t. The background change is controlled by the weight η μ , which determines how much the background of the next frame should differ from the present one. Calculating the binary difference between the current frame F t (i, j) and the estimated background E t (i, j) yields the foreground by If the difference between each pixel coordinate of F t (i, j) − E t (i, j) is smaller than the predefined threshold value η α , we will consider that pixel point to be zero (dark) in image K(i, j). The threshold η γ may be automatically determined by where n h , n v are the number of image pixels in horizontal and vertical directions of the image and η λ is the coefficient of inhibition. The value of η α ranges between 10 and 30. The images may still have some noise. To get rid of the unwanted elements "Gaussian blurring" and "morphological closure" are applied by opencv image processing toolbox where Gaussian kernel acts as a low pass filter to eliminate high frequency components and morphological closure helps to fill small gaps in the images. After the stimulus and insect have been segmented, we calculate the area size, two-dimensional center position (centroid), shapes, and other metrics. Occasional frames may be missed owing to overlaps or occlusion. Using the discrete Kalman filter [43], each insect centroid position can be predicted when there are any missing frames. Detailed information about the Kalman filtering method is described in appendix (section Kinematic filtering). These 2D centroids are required in the next stage. 2D to 3D conversion. 3D position from the 2D measured positions of all cameras can be calculated by camera calibration matrices. The intrinsic and extrinsic calibration parameters for each of camera are found using the Svoboda multi camera calibration toolbox [76]. During the Matlab execution of this toolbox, a red laser light in the working volume must be moved. Using these calibration matrices and previously tracked 2D points, we employ the linear triangulation approach to obtain 3D points. For each camera model j, the linear triangulation equation is stated as follows v j i P j 3 À P j 2 u j i P j 3 À P j 1 :::: :::: where ðu j i ; v j i Þ is 2D point i of j th camera, O i is the 3D point and P j n is the n th row of the calibration matrix P 3×4 of camera j. Combining perspective projection expressions for multiple cameras yields a homogeneous system of linear equations; at least two cameras are required to solve this system of equations via singular value decomposition. The combination of 2D points acquired from the cameras can be used to calculate 3D points. The 3D point then is re-projected back to the 2D point by the calibration matrices. A tolerance η ψ on the re-projection error between the true and re-projected 2D points is used to specify the desired precision.

System identification
In the frequency domain, tracking trajectories may be described by gain, phase, and coherence. The gain depicts the insect's amplitude in relation to the stimulus in each frequency component. To translate the trajectory from time domain to frequency domain, we apply the Fourier transform (FFT). Tracking error e(s) may be computed by reference to ideal tracking (gain 1, phase 0) as where G(s) is the transfer function containing gain and phase for each frequency s = jω point. Coherence γ 2 (s) is calculated from the spectral power density of an input and output signal pair (x, y) by the following equation as where S xy (s) denotes the cross spectral power density and S x (s) and S y (s) represent the auto power spectral density of the stimulus (input) and bee (output) coordinates, respectively. The coherence between the stimulus and the bee trajectory is used to determine the degree of linear connection throughout the frequency range with 1 denoting a purely linear system relationship and 0 denoting no linear relationship. For this investigation, we assumed a linear relationship between the stimulus and the bee trajectories up to a coherent frequency with a coherence value greater than 0.7. The FFT transform outputs a discrete time frequency domain signal up to the Nyquist frequency ( 1 2dt ). To improve the resolution over a targeted frequency range identified from the coherence, Chirp Z transform (CZT) [77] was applied. CZT is a generalization of discrete Fourier transform (DFT) and samples the Z plane via spiral arcs (straight lines in the S plane), whereas the DFT samples the complex Z plane at equally spaced points along the unit circle. We consider this CZT frequency domain data of the stimulus D 1 (s) and bee D 2 (s) trajectories to examine the flight dynamics of the tracking trajectories. We conduct the system identifica- where y is the true output, y m is the model's predicted output and n is number of data points.
To find the FPE, we have considered an ARX (auto-regressive with exogenous input) model as yðtÞ þ a 1 yðt À 1Þ þ :::: þ a n yðt À n a Þ ¼ b 1 uðt À 1Þ þ ::: þ b n uðt À n b Þ::: þ wðtÞ ; ð11Þ where n a , n b , w(t) are the number of poles, number of zeros and white noise term respectively. The estimated parameters are θ = [a 1 , a 2 , ...., a n b 1 , b 2 , ...., b n ] T and regression matrix is where L is the number of values in the estimate data set and the number of estimated parameters is represented by d.

Swarm model
Our experiment extracts single bee tracking dynamics, and we build a hypothetical swarm of visually interconnected insects from those delays and transfer functions. Each swarm agent individual has a distinct communication latency, consistent with the individual reaction time measurements. A delay-based factor among the insects may impact group behaviors, and we included a coupling strength to allow one to consider multiple swarm shapes.

Swarm dynamics
We consider a second order (position and velocity) dynamic system model for N agents communicating with each other agents with some heterogeneous processing delays. These delays are heterogeneous i.e, each agent has distinct delay. The dynamic model based on 2D or 3D position x i and velocity v i may be written in either as where r i A a (x ij , τ ij ) and r i A r (x ij , B r , C r ) are an attractive and repulsive potential respectively, specified as Here, τ ij is the delay from agent i to agent j. x ij is inter-agent distance of agent i to agent j with some time delay τ ij and is denoted by . ρ is inter-agent coupling strength and we assume each agent goes to r distance from the origin. B r and C r are the amplitude and applied distance of repulsive potential. β is the coefficient of friction. For notational convenience, let X i ¼ P N i¼1;i6 ¼j ðx i ðtÞ À x j ðt À t ij ÞÞ. The magnitude of the movement I is written as and the attractive potential gradient is where cosðyÞ ¼ X i ð1Þ The center of mass of the swarm R C can be defined as To find the position stability of the swarm, the swarm center norm can be taken as kR C k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi , where R CX , R CY are the X and Y coordinates of the swarm center respectively.
For the simulation described in Analysis & simulation, each agent's position output x i then passes through its identified transfer function G ei (s). If the identified transfer function is G ei ¼ a n s 2 þa 1 sþ::þa 0 1þb 1 sþ:::b n s n , then the filtered position � x i ðtÞ will be, � x i ðtÞ ¼ a 0 x i ðtÞ þ ::: þ a n x i ðt À nÞ À b 1 � x i ðt À 1Þ À :::: À b n � x i ðt À nÞ: ð22Þ Linear stability. Stability at long range is characterized by the attraction potention (the underlying mathematical theory to show this is included in Appendix section "Long range attraction"). Linear analysis characterises the system's local stability properties. The long range dynamics (attraction potential) may be rewritten as d dt and v 0 = 0, A a (x 0 ) = 0 at the equilibrium point (x 0 , v 0 ). This occurs when X i = r if ρ < 0 and X i ¼ r; r � ffi ffi ffi r p if ρ > 0. The system eigenvalues may be obtained via the Jacobian matrix J, given by where J 21 can take on three values depending on ρ and equilibrium points X i . The constant r in this case represents the distance from the origin.
When J 21 = ρ, the eigenvalues will be l ¼ À b=2 � 1 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi b 2 þ 4r p and ρ < 0 is sufficient to ensure the eigenvalues have negative real part and the equilibrium is stable. When ρ > 0, at least one positive real eigenvalue exists and the equilibrium is unstable. For J 21 = −2ρ, the equilibrium is ffi ffi ffi r p and the eigenvalues are l ¼ À b=2 � 1 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi b 2 À 8r p . Stability of this equilibrium requires β > 0 and ρ > 0.
The attraction potential can be used as a pitch fork bifurcation equation with ρ as the bifurcation parameter and forms different stable conditions depending on the sign of ρ. This is a supercritical pitchfork bifurcation in which ρ bifurcates into two local equilibria from a single equilibrium. Table 1 shows the linear stability based on the sign of ρ. This stability analysis finds that there are three possible swarm patterns that may be generated using the bifurcation parameter. The theoretical analysis does not consider the effects of delays. Our expectation is that the time delays as a whole stabilize or destabilize the formation of the swarm shape which is shown in the simulation.

Results and discussion
An example insect trajectory and its analysis is shown in the experimental section. We have done similar approach to perform system identifications in all other insects. The simulation section demonstrates how changing the bifurcation parameter may produce various swarm configurations. Later, the simulation result demonstrates the impact of identified time delays and transfer functions on the swarm center position.

Experimental
The input stimulus was generated by eight different sinusoidal signals whose frequency range is 0.1-1.7 Hz and a bee's response is shown in Fig 4. Using the tracking system 3D position data was found. The parameters of the tracking system used in the experiment are shown in Table 2. The bee tracked the light stimulus intermittently.    The bee's delay in responding to the target is indicated by a negative phase (lag).Using coherence we can identify tracking performance over a frequency range. The coherence plot of Fig 6a shows that the input and output have a strong linear relationship up to 0.1 Hz (as quantified by γ 2 > 0.7) and its error e(s) is small (<5dB) in this range. The number of points in this range, however, is less than 20. To improve the resolution in the region of high coherence, the CZT transform was used. The CZT frequency upper limit, shown in Fig 6b, was set to 0.2 hz to increase resolution in the dominant frequency range. The frequency response function H CZT ðsÞ from the CZT was then used to find the best matched transfer function model. As seen by the FIT error statistics presented in Table 3, a 2 pole, 1 zero transfer function with 21 ms processing latency (transport delay) was the best fit for this example trajectory. The  When a transfer function model was fit to the frequency range of the experimental frequency responses having high coherence, best-fit model in terms of pure tracking delay and linear approximation was seen in Table 4. When identifying the system dynamics, the uncertainty of

Analysis & simulation
When coupling strength ρ is used as a pitch fork bifurcation parameter, a visually-interacting swarm is able to produce both stable and unstable modes, as quantified by linear stability analysis. The achievable parameter values and pattern shapes are summarized in Table 5, which shows the effect of varying ρ on resulting formations (the variation in number of agents N was for visualization and did not affect on the structure). The visually-interconnected swarm simulations showed a diversity of achievable patterns consistent with theoretical predictions in which ρ was varied to establish the cluster form. The simulations may be conducted in 2D and 3D, with no theoretical impact. Here, we present 2D simulations for computational and presentation simplicity. Parallel 3D cluster formation is illustrated in the appendix (e.g., Fig 8). In both cases, the existence of stable behavior depends on coupling strength ρ. The choice of ρ = −4 creates a single ring structure. As ρ decreases to 1.5, the ring's stable equilibrium becomes unstable and it bifurcates into two rings. Simulations conducted with the same initial positions

PLOS ONE
illustrate this effect, as shown for three cases in S2 Video and (Fig 9b-9d). According to the Eq (20), the simulated agents' position and velocity meet the consensus conditions after a period of time, as seen in (Fig 9e and 9f).

Effect of identified delays
We investigated the effect of delays on the swarm center. The more naturalistic cluster shape was used to evaluate the swarm's center of mass (barycenter) position stability under varying interaction delays. The simulation was run under different positive interaction delays based on a Gaussian distribution N ðt > 0jm; sÞ, which is defined as N ðt > 0jm; sÞ ¼ 1 where mean μ and standard deviation σ were varied. Plotting the mean μ and standard deviation σ illustrates two distinct stable and unstable areas. As illustrated in Fig 10a, the black points represent the unstable behavior and the green points are for the stable behavior. A third order polynomial was sufficient to describe the boundary between the stable and unstable regions. In this simulation, interaction delays for each participant in a visually interconnected swarm is assigned a delay randomly from the measured distribution. Because the measured delay histogram approximates a normalized Gaussian distribution with mean μ = 22 and standard deviation σ = 40 ms located in the unstable region of Fig 10a, our analysis indicates the simulation using measured delays will have an unstable center of mass. We also wanted to simulate using the actual experimental delays. Fig 10b shows the simulation result using the experimental measured delays for 10 seconds. In Fig 10c, the norm of swarm centers ||R c || are shown for all three cases: no delay, Gaussian delays in stable region, and experimental delays. The ||R c || for measured delays shows position instability, while the delay free case and delays within stable area of Fig 10a show position stability of ||R c ||.

Effect of both identified delays and transfer functions
We also conducted a simulation including both the identified individual transfer functions as well as the delays, which shows cluster shape formation in the stable region shown in Fig 11. The inclusion of individual agent transfer functions changes the location of the swarm center rather than the overall shape of the swarm. Finally we applied the identified transfer functions and delay variation and its swarm center position is shown in Fig 10c.

Conclusion
In this study, individual honeybees tracked a visual light stimulus and the visuomotor delay in their closed loop tracking was measured. We developed a real-time camera-based tracking system called VISIONS that track honeybee 3D position and induced them to follow a moving light target. A system identification technique was applied to identify the closed loop tracking dynamics between light stimulus motion and insect body motion, and quantify the delay between the stimulus and animal trajectories, separating the effects of open loop plant (locomotion) from visuomotor feedback dynamics. The measured honeybee sensorimotor delays were used to find a delay distribution across population, showing that insect visual sensorimotor delays in a tracking task are heterogeneous across population.
To understand the implications of the measured delays on visual communication and identified dynamic systems, we then integrated the measured delays and dynamic systems into a visually interacting swarm model. Analysis on this model indicates the range of achievable swarm patterns (cluster, ring, double-ring) and conditions needed for center of mass's position stability of each mode, and simulation illustrates these achievable behaviors and stability regions for both theoretical and measured delays. The analysis and simulation indicate that while the processing delays measured in solitary conditions support three relative shapes, these delays lie in a region associated with an unstable center of mass position and are thus sufficient to support coordinated relative motion but not center of mass position stabilization. This finding suggests that visually interconnected peer insects may be able to to achieve relative configurations but that the measured visual interconnection structure does not support stabilizing the swarm's overall center of mass position. An important distinction is that this study considers conspecific peers and not the effect of other visual targets or stimulus, which may play a role in the group's position. In the absence of external targets or a form of delay modulation (compensation), this analysis and lack of barycenter position stiffness suggests that swarm center of mass position may drift, an important result to completing swarm theory descriptions and informing experiments that investigate solitary and swarm motions in flying insects.
The prediction matrix P − (t) and integratedx À ðtÞ are applied as In the correction step, the Kalman gain K is obtained as Finally, we update the P + (t) andx þ ðtÞ by x þ k ¼x À k þ K k ðy k À Cx À k Þ : ð39Þ

Long range attraction
At long distance the repulsive potential have negligible effect on the swarm model. The distance between every agent i to other agent j is taken by U = |x i (t) − x j (t − τ ij )|. From Eq (14) Thus at long distance the repulsive force vanishes, and the attraction potential characterizes long range stability analysis. 3D simulation. 3D Cartesian coordinates in Eq 18 must be considered in order to simulate the swarm model in three dimensions. Here, a 3D simulation of 50 agents is illustrated using the parameters ρ = −4, r = 0, B r = 1, C r = .5, β = 10, and random initial conditions. This cluster swarm formation shows that the swarm model works in both 2D and 3D environments.